%%%%%%%%%%%%%%%
% Housekeeping
%%%%%%%%%%%%%%%

clc;
clear                                   all;
close                                   all;
beep                                    off;
rng.Seed                                = 110887;
set(groot,'defaultAxesFontName', 'Garamond')
set(groot,'defaultTextFontName', 'Garamond')
set(groot,'defaultAxesTickLabelInterpreter','latex'); 
set(groot,'defaultLegendInterpreter','latex');
cd(pwd); 
addpath(strcat(pwd,'/auxiliary_functions/'));
project_folder                          = fileparts(fileparts(pwd));
root_folder                             = fileparts(project_folder );

%%%%%%%%%%%%%%%
% Model setup
%%%%%%%%%%%%%%%

% Structural parameters
pars.r                                  = -1/12*log(0.70);                                                              % Monthly model with annual discount rate of 0.70 
pars.k                                  = 0.163;                                                                        % Upper bound on adjustment intensity
pars.M_c                                = 1;                                                                            % Cash stock average
pars.M_e                                = 0.97;                                                                         % Electronic money    
pars.C                                  = 0.060;                                                                        % Strength of complementarities           
pars.mean_reversion                     = 120;                                                                          % Shock persistence, expressed as the number of days to get back to within 10% of initial value
pars.theta                              = -log(1-0.90)/( pars.mean_reversion/30 );                                      % Shock persistence, expessed as OU parameter
%pars.theta                              = 0.2;
pars.sigma                              = 0.060;                                                                        % Standard deviation of money shock    
pars.T                                  = 12*10;                                                                        % Date after which mean-reversion goes to zero
pars.tIRF                               = 12*2;                                                                         % IRF horizon

% Discretization parameters
Nu.NX                                   = 2e1+1;                                                                        % # of gridpoints for X
Nu.N                                    = 2e1;                                                                          % (1/2)*( # of gridpoints for M - 1 )
Nu.Nmin                                 = ceil(pars.theta/pars.sigma^2* ...
                                               ((pars.r+pars.k+pars.theta)/(pars.r+pars.k))^2* ...
                                              max((pars.M_c-pars.M_e)^2,(pars.M_e + pars.C - pars.M_c)^2)); 
Nu.N                                    = max(Nu.N,Nu.Nmin); 
Nu.Deltat                               = 1/(Nu.N*pars.theta);                                                          % Time period, expressed as fraction of a month
Nu.theta_target                         = pars.theta;
Nu.solutionmode                         = 1;

%%%%%%%%%%%%%%%%%%%%%%
% Construct thresholds
%%%%%%%%%%%%%%%%%%%%%%

S                                       = 0.2;
NX                                      = 1e2;
X                                       = linspace(0,1,NX)';
X_underline                             = NaN(size(X));
C_overline                              = NaN(size(X));
X_overline                              = NaN(size(X));
C_underline                             = NaN(size(X));
a_LB                                    = (1+pars.theta/(pars.r+pars.k))*( (pars.M_c-pars.M_e)/(pars.M_c) );
options                                 = optimset('Display','off');
for n=1:NX
    Xn                                  = X(n);
    if Xn < 1
       RHS_LB                           = (1-Xn).^(pars.theta/pars.k)*a_LB*pars.k/S;
       f_LB                             = @(X) pars.theta*(1-X).^(pars.theta/pars.k-1) - (pars.theta-pars.k)*((1-X).^(pars.theta/pars.k))-RHS_LB; 
       X_underline(n)                   = fzero(f_LB,[Xn 1]);
       b_LB                             = S*pars.theta/pars.k*(1-X_underline(n)).^(pars.theta/pars.k-1).*(1-Xn).^(-pars.theta/pars.k);
       C_overline(n)                    = (pars.r+2*pars.k)/(pars.r+pars.theta+pars.k)*pars.M_c*b_LB;
       RHS_UB                           = RHS_LB;
       f_UB                             = @(X) (1+pars.k/(pars.r+pars.k))*pars.theta*(1-X).^(pars.theta/pars.k-1) - (pars.theta-pars.k)*((1-X).^(pars.theta/pars.k))-RHS_UB; 
       X_overline(n)                    = fzero(f_UB,[Xn 1]);
       b_UB                             = S*pars.theta/pars.k*(1-X_overline(n)).^(pars.theta/pars.k-1).*(1-Xn).^(-pars.theta/pars.k);
       C_underline(n)                   = (pars.r+2*pars.k)/(pars.r+pars.theta+pars.k)*pars.M_c*b_UB;
    else
       X_underline(n)                   = 1;
       C_overline(n)                    = (pars.r+2*pars.k)/(pars.r+pars.k)*(pars.M_c-pars.M_e);
       X_overline(n)                    = 1;
       C_underline(n)                   = (pars.M_c-pars.M_e);
    end
end

oC                                      = C_overline;
uC                                      = C_underline;

%%%%%%%%%%
% Plots
%%%%%%%%%%

fignum                                  = 1;
color1                                  = [0, 0.4470, 0.7410];
color2                                  = [0.6350, 0.0780, 0.1840];
gray                                    = [0,0,0]+0.3;
gray2                                   = [0,0,0]+0.5;
gray3                                   = [0,0,0]+0.7;
gray4                                   = [0,0,0]+0.95;
color3                                  = [0.9290, 0.6940, 0.1250];
color4                                  = [0.4660, 0.6740, 0.1880];
linewidth                               = 2;

% theta <= k
ff3 = figure(fignum);
fignum = fignum+1;

    oC = oC(end)*ones(size(oC));
    uC = uC(end)*ones(size(uC));
    p1=plot(X,uC,'-.','LineWidth',linewidth,'Color',gray2);
    hold on
    p2=plot(X,oC,'--','LineWidth',linewidth,'Color',gray2);
    hold on;
    ax = gca;
    ylim([0 0.15])
    rna1 = fill([X' flip(X)'],[(oC') (flip(uC)')],gray3,'LineStyle','none');
    rna1.FaceAlpha = 0.5;
    rna2 = fill([X' flip(X)'],[(oC') (ax.YLim(2)*ones(size(uC)))'],color1,'LineStyle','none');
    rna2.FaceAlpha = 0.15;
    rna3 = fill([X' flip(X)'],[(ax.YLim(1)*ones(size(uC)))' (flip(uC)')],color2,'LineStyle','none');
    rna3.FaceAlpha = 0.15;
    ylabel({'$C$ \quad'},'interpreter','latex','Rotation',0);
    xlabel('$X_0$','interpreter','latex');
    xlim([0 1]);
    xticks([0 1]);
    set(gca,'ytick',[]);
    set(gca,'yticklabel',[]);
    set(findall(gcf,'-property','FontSize'),'FontSize',20); 
    text(0.5,0.5*uC(ceil(NX/2)),{'$\hat{t}(X_0) < +\infty$'},'color',color2,'interpreter','latex','HorizontalAlignment','center','fontsize',15);  
    text(0.5,uC(ceil(NX/2))+0.4*(oC(ceil(NX/2))-uC(ceil(NX/2))),{'$\hat{t}(X_0) \leq +\infty$'},'color',gray2,'interpreter','latex','HorizontalAlignment','center','fontsize',15);  
    text(0.5,2*oC(ceil(NX/2)),{'$\hat{t}(X_0) = +\infty$'},'color',color1,'interpreter','latex','HorizontalAlignment','center','fontsize',15);  
    legend([p2,p1],{'$\overline{C}(X_0)$','$\underline{C}(X_0)$'},'interpreter','latex','location','Northeast','fontsize',15);

    % Print
    set(ff3,'Units','inches');
    screenposition                              = get(ff3,'Position');
    set(ff3,'PaperPosition',[0 0 screenposition(3:4)],'PaperSize',[screenposition(3:4)]);
    ff3.PaperPositionMode                       ='auto';
    set(ff3,'position',[0,0,8,7])
    exportgraphics(ff3,[root_folder '/output/figures/' 'Figure_H_16.pdf'],'BackgroundColor','none')